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1  Introduction 


Hearing  impairment  remains  the  primary  disability  among  military  person¬ 
nel.  Sound  pressure  levels  caused  by  proximity  to  aircraft  engines  or  impulse 
noise  from  large  caliber  weapons  may  easily  exceed  the  pain  threshold  value 
of  100  dB. 

The  focus  of  this  project  is  to  to  develop  a  reliable  numerical  model  for 
investigating  the  bone-conducted  sound  in  the  human  head.  The  problem 
is  difficult  because  of  a  lack  of  fundamental  knowledge  regai’diiig  the  trans¬ 
mission  of  acoustic  energy  through  non-airborne  pathways  to  the  cochlea.  A 
fully  coupled  model  based  on  the  acoustic/clastic  interaction  problem  with 
a  detailed  resolution  of  the  cochlea  regton  and  its  interface  with  the  skull 
and  the  air  pathways,  should  provide  an  insight  into  this  fundamental,  long 
standing  research  problem. 

The  project  builds  on  an  interaction  of  experts  in  numerical  wave  prop¬ 
agation  -  Drs.  Elizabeth  and  Marek  Bleszynski  from  Monopolc  Research 
with  a  team  at  the  University  of  Texas  headed  by  Dr.  Leszek  Demkowicz 
and  including  two  experts  on  wave  propagation  and  hearing  science:  Dr. 
Mark  Hamilton  and  Dr.  Greg  Champlin. 

2  The  Head  Problem 

In  this  section  we  review  shortly  specifics  of  the  head  problem.  The  prob¬ 
lem  falls  into  the  category  of  general  coupled  elasticity /acoustics  problems 
discussed  in  Appendix  A  with  a  few  minor  modifications.  The  domain  in 
which  the  problem  is  defined  is  the  interior  of  a  ball  including  a  model  of 
the  human  head,  and  it  is  split  into  an  acoustic  part  Qa^  and  an  elastic  part 
De.  Depending  upon  a  particular  example  we  shall  study,  the  acoustic  part 

includes: 

•  air  surrounding  the  human  head,  bounded  by  the  head  surface  and  a 
truncating  sphere;  this  part  of  the  domain  may  include  portions  of  air 
ducts  leading  to  the  middle  car  through  mouth  and  nose  openings; 

•  cochlea, 

•  an  additional  layer  of  air  bounded  by  the  truncating  sphere  and  the 
outer  sphere  terminating  the  computational  domain,  where  the  equa¬ 
tions  of  acoustics  arc  replaced  with  the  corresponding  Perfectly  Matched 
Layer  (PML)  modification. 
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The  elastic  part  of  the  domain  hiciudes: 

•  skull, 

•  tissue. 


By  the  term  tissue  we  understand  here  all  parts  of  the  head  that  arc  not 
occupied  by  the  skull  (bone)  and  the  cochlea.  This  includes  the  thin  layer  of 
the  skill  and  the  entire  interior  of  the  head  with  the  brain.  We  will  assume 
that  the  clastic  constants  for  the  whole  tissue  domain  are  the  same.  As  the 
viscosity  constant  for  the  tissue  is  four  orders  of  magnitude  smaller  than 
that  for  the  skull,  an  alternative  approach  would  be  to  model  the  tissue  as 
an  acoustical  fluid,  neglecting  the  shear  waves  in  the  brain. 

The  acoustic  wave  is  reprcsctitcd  as  the  sum  of  an  incident  wave 
and  a  scattered  wave  p.  Only  the  scattered  wave  is  assumed  to  satisfy  the 
radiation  ( Sommer feld)  condition, 


^  +  ikp  e  L^(R^) 
or 


(2.1) 


The  different  types  of  boundaries  discussed  in  Appendix  A  reduce  only  to 
the  interface  between  the  clastic  and  acoustic  subdomains,  and  the  outer 
Dirichlct  boundary  for  the  acoustic  domain.  Material  interfaces  between 
the  skull  and  tissue,  as  well  between  the  air  and  the  PML  air  do  not  require 
any  special  treatment. 

The  final  formulation  of  the  problem  has  the  form  A. 11,  with  the  bilinear 
and  linear  forms  defined  as  follows. 


bee{u,  v)  = 

f»ae(p,  W)  = 


?)  — 


baaip,  9)  = 


Uv)  = 
UiQ)  = 


/  {EijklUk,tVij  -  Pau'^tliVi)  dx 
Jiit 

/  pVn  dS 
Jrt 

-iJpj  /  dS 

Jr, 

I  {Vp'Vq  —  k^pq)  dx 
Jna 

-  f 

Jr, 


(2.2) 
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material 

£;[MPa] 

1/ 

p[kg/m^ 

r.p[m/s] 

Ca[m/s] 

tissue  (brain) 

0.67 

0.48 

1040 

75 

15 

skull  (bone) 

6500 

0.22 

1412 

2293 

1374 

cochlea  (water) 

1000 

1500 

air 

1.2 

344 

Table  1:  Material  constants  and  speed  of  sound  for  comprcssional  and  shear 
waves. 


A  symmetric  formulation  is  enabled  by  dividing  the  equations  of  acoustics 
with  factor 


u)  = 
baeiP,  V)  = 
beaiu,q)  ~ 
f^aa{Pt  q)  ~ 

Uiv)  = 

Uq)  = 


/  {EijkiUk,iVi,j  -  p^ui'^UiVi)  (lx 

L 


pUn  dS 


-  /  UnqdS 

hi 

-  /  dS 

h, 


(2,3) 


Notice  that  we  refer  to  the  outer  normal  unit  vector  n  always  locally,  i.c. 
in  the  formula  for  the  coupling  bilinear  form  60^  involving  elasticity  test 
functions  u,  versor  n  points  outside  of  the  clastic  domain ,  whereas  in  the 
formula  for  the  coupling  bilinear  form  involving  acoustic  test  functions 
q,  versor  n  points  outside  of  the  acoustic  domain*  The  normal  components 
Vn  and  Utj  present  in  the  coupling  terms  arc  thus  opposite  to  each  other,  and 
the  formulation  is  indeed  symmetric. 


2.1  Material  and  temporar  scales 

The  material  data  are  summarized  in  Tabic  1.  For  clastic  materials,  the 
speed  of  comprcssional  waves  and  shear  waves  is  given  by  the  formulas: 


2  E{1  “  £/)  2  _  ^ 

""p  ^  (i+i/){i-2i/)p’  *  “  2(r+vj“p 


(2.4) 
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It  is  illuminating  to  begin  with  an  estimate  of  the  magnitude  of  all  terms 
entering  the  equations.  There  are  three  scales  to  the  problem. 

Element  scale:  h  representing  a  typical  element  length.  Assume  h  ^  l[em] 
-  10-2  [m]. 

Elastic  displacement  in  skull  scale:  b  representing  the  expected  magni¬ 
tude  of  the  acoustical  displacements  within  the  skulL  Assuming  the  air 
pressure  at  the  level  of  lOOfPa  =  N/rn^j  (threshold  of  pain),  a  simple 
model  of  a  pressure  loaded  elastic  shell  of  thickness  1-2  [cm]  responds 
with  a  displacement  of  order  b  =  0,l[mmj=  10“^  [m]. 

Pressure  scale:  po  to  be  assumed. 

The  following  table  represents  the  order  of  magnitude  of  all  terms  entering 
the  formulation  for  the  coupled  problem  involving  the  skull  and  the  air  only, 
for  an  average  frequency  of  /  =  200 [Hz]. 

elastic  stiffness  term  = 
elastic  mass  term  = 
coupling  term  ^ 

acoustic  stiffness  term  = 
acoustic  mass  term  = 


fib^h  ^  1 


^  10  ^ 

PQbh^  10”®po 


1 


=  lo-'VS 


(2.5) 


Assuming  the  pressure  factor  in  the  range  of  po  ^  10®,  we  observe  that  all 
but  one  terms  in  the  formulation  vary  by  two- three  orders  of  magnitudes 
only.  Notice  that  the  assumed  scales  for  displacement  and  pressure  have 
actually  nothing  to  do  with  the  expected  displacement  and  pressure  levels. 
Assuming  that  displacements  are  of  order  1  would  yield  the  same  result. 
For  convenience,  wc  have  also  assumed  the  same  order  of  magnitude  for  test 
functions.  Again,  due  to  the  linearity  of  the  equations,  the  choice  docs  not 
matter.  The  five  orders  of  magnitude  difference  for  the  elastic  mass  term 
(the  subwavciength  regime)  should  not  present  a  problem  in  double  precision 
computations. 

In  conclusion,  the  advantage  of  the  coupled  formulation  ^'mixing”  pri¬ 
mary  variables  {clastic  displacement)  with  dual  variables  (pressure),  is  the 
possibility  of  adjusting  the  scaling  to  yield  terms  of  the  same  order.  The 
scaling  arguments  have  been  confirmed  by  monitoring  pivots  reported  by 
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frontal  solver  when  solving  examples  discussed  in  the  next  section.  Indeed, 
the  pivots  varied  by  the  three-four  orders  of  magnitude  only. 


2,2  PML  modification 


In  the  PML  part  of  the  acoustical  domain,  the  bilinear  form  is  modified 
as  follows: 


Ka{p,q)= 


dpdg  dp  dq 


dpdq\ 


dr  dr  dip  dip 


sin^  ip  dd  dd  ) 


sin  ipdrdipdO 


(2.6) 

Here  r,  0  denote  the  standard  spherical  coordinates  and  z  =  z{r)  is  the 
PML  stretching  factor  defined  as  follows, 


z{r)  =  {1 


i 

k 


{2.7) 


Here  a  is  the  radius  of  the  truncating  sphere,  h  is  the  external  radius  of 
the  computational  domain  {6  ^  a  is  thus  the  thickness  of  the  PML  layer),  i 
denotes  the  imaginary  unit,  k  is  the  acoustical  wave  number,  and  r  is  the 
radial  coordinate.  In  computations,  all  derivatives  with  respect  to  spherical 
coordinates  arc  expressed  in  terms  of  the  standard  derivatives  with  respect 
to  Cartesian  coordinates.  In  all  reported  computations,  parameter  q  =  5. 
For  a  detailed  discussion  on  derivation  of  PML  modifications  and  effects  of 
higher  order  discretizations  see  [9]. 


3  Numerical  Examples 

We  refer  to  Appendix  B  for  details  on  our  finite  elcnient  technology  and  to 
Appendix  C  for  a  description  of  the  parallel  linear  solver  used  in  Phase  I  of 
this  project. 


3,1  Verification 

The  code  has  been  verified  by  implementing  so-called  manufactured  solu¬ 
tions.  This  is  a  standard  technique  in  finite  elements.  We  assume  an  ana¬ 
lytical  solution  of  any  form  (the  manufactured  solution),  and  use  the  differ¬ 
ential  equations  for  both  acoustics  and  elasticity  parts,  boundary  and  inter* 
face  conditions,  to  compute  the  corresponding  volume  forces  and  boundary 
fluxes.  The  verification  is  invaluable.  By  assuming  a  solution  that  can  be 


5 


Figure  1:  Verification  of  the  code  using  a  manufactured  solution-  Pressure 
distribution  along  a  section  passing  through  the  origin  and  parallel  to  the 
x^axis  for  (a)  unit  sphere  testj  and  (b)  Example  1  with  unit  material  data. 
Numerical  and  exact  (manufactured)  solutions  are  indistinguishable. 


reproduced  exactly  with  the  FE  shape  functions,  we  do  know  that  the  cor¬ 
responding  error  must  be  equal  to  machine  zero,  in  our  case  values  around 
10“^"^.  Any  values  bigger  than  these,  indicate  a  bug  in  the  code.  Values 
around  10"^  indicate  a  loss  of  double  precision.  In  our  casc^  due  to  the  use 
of  linear  elements  in  the  acoustic  domain  and  linear  or  quadratic  elements 
in  the  elastic  pai't,  we  can  assume  any  linear  variation  for  pressure,  and 
any  linear  (quadratic)  variation  for  displacement  vector.  To  verify  the  codc^ 
we  have  used  a  simple  example  of  domain  consisting  of  a  unit  sphere,  sur¬ 
rounded  with  a  unit  layer  of  air,  and  a  PML  layer  of  thickness  equal  to  two 
units-  The  sphere  was  meshed  with  just  eight  octant  tetrahedra,  and  the 
air  layer  with  8x3  tetrahedra  obtained  by  splitting  eight  prisms,  each  into 
three  tetrahedra.  The  PML  layer  was  modeled  with  two  layers  of  prismatic 
elements  with  arbitrary  order  p  <  4  in  the  radial  direction.  All  material 
data  were  set  to  0(1)  values.  A  typical  result  of  the  verification  for  a  mesh 
with  quadratic  elements  in  the  clastic  domain,  is  shown  in  Fig,  L  The  same 
verification  technique  was  then  use  to  verify  each  data  set  and/or  the  three 
different  solvers  used  for  the  project:  frontal  solver,  MUMPS,  and  the  paral¬ 
lel  nested  dissections  multifrontal  solver.  For  small  data  sets,  an  additional 
verification  is  done  by  comparing  results  obtained  with  the  different  solvers. 

3-2  Phase  I  Examples 

Two  numerical  cxainpics  have  been  considered  in  the  first  phase  of  the 
project. 
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3*2,1  Scattering  of  a  plane  wave  on  an  elastic,  multilayer  sphere 

In  this  model,  the  domain  consists  of  four  concentric  spheres.  The  most 
inner  sphere  is  filled  with  an  clastic  material  with  data  corresponding  to 
human  brain.  The  first  layer  is  also  clastic  with  constants  corresponding 
to  human  skull.  The  second  layer  corresponds  to  air,  and  the  last  one  to 
the  PML  air.  The  incident  wave  is  assumed  in  the  form  of  a  plane  wave 
impinging  from  the  top, 

p'"'"  =  e  =  (0,0, -1),  po  =  l[Po]  (3.1) 

The  test  problem  is  being  solved  with  frequency  /  =  200  Hz.  The  precise 
geometry  data  are  as  follows: 

brain  r  <0.1  m 

skull  0.1  m  <  r  <.125  m 

air  0,125  m  <  r  <  0.2  m 

PML  air  0.2  m  <  r  <  0.3  m 


Three  tetrahedral  meshes  for  the  interior  tissue  ball,  of  radius  lOon  ,werc 
generated  using  the  simple  MATLAB  code  distmesh,  described  in  [6].  The 
surface  of  this  mesh  was  then  manually  extended  to  generate  a  skull  annulus 
of  thickness  2. Son,  surrounding  air  of  thickness  7.5cm,  and  a  PML  of  thick¬ 
ness  10cm.  The  problem  was  solved  on  three  meshes  shown  in  Figures  2,  3, 
and  4.  ,  We  will  refer  to  them  as  ''smair\  '‘big”,  and  ‘‘huge”  meshes.  For 
all  runs  discussed  for  this  example,  we  have  used  the  MUMPS  solver. 

Fig.  5  displays  the  distribution  of  the  real  part  of  the  pressure  over  plane 
y  =  0  passing  through  the  origin.  It  looks  “good”.  Unfortunately,  a  similar 
picture  for  the  imaginary  part  reveals  a  severe  instability  in  the  “tissue” 
region.  To  double  check  the  VTK  graphics,  we  have  displayed  the  results 
across  a  vertical  section  passing  through  the  origin.  The  results  are  shown 
in  Fig.  6,  In  order  to  access  the  problem,  we  have  run  the  same  example 
but  with  the  Young  modulus  for  the  tissue  domain  increased  by  two  orders 
of  magnitude,  i,e.  E  =  67  [MPa].  The  corresponding  results  are  shown  in 
Fig.  6,  Figure  8  displays  the  same  pressure  in  dB. 

Figures  9  and  10  display  the  same  pressure  but  this  time  obtained  on 
the  big  and  huge  meshes.  The  values  arc  at  the  same  level  which  indicates 
a  converged  solution. 

Finally,  Figures  11  and  12  display  the  distribution  of  real  and  imaginary 
parts  of  the  pressure  over  the  j/  =  0  section  obtained  on  the  big  mesh. 
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Figure  2;  The  “small”  mesh  used  to  solve  the  multilayer  sphere  problem. 
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Figure  3:  The  “big”  mesh  used  to  solve  the  multilayer  sphere  problem. 
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Figure  4:  The  “huge”  mesh  used  to  solve  the  multilayer  sphere  problem. 
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Figure  5:  The  concentric  spheres  problem.  Small  mesh.  Plot  of  real  part  of 
pressure  on  plane  y  —  0  passing  through  origin. 


Figure  6:  The  concentric  spheres  problem.  Small  mesh.  Plots  of  real  and 
imaginary  part  of  pressure  along  the  vertical  seetion. 
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Figure  7:  The  concentric  spheres  problem.  Small  mesh.  Plots  of  real  and 
imaginary  part  of  pressure  along  the  vertical  section  for  the  case  of  a  ^^stifFer” 
tissue. 
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Figure  8:  The  concentric  spheres  problem.  Small  mesh.  Plot  of  pressure 
along  the  vertical  section  in  dB  for  the  case  of  a  *'stiff'’  tissue. 
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Figure  9:  The  concentric  spheres  problem.  Big  mesh.  Plot  of  pressure  along 
the  vertical  section  in  dB  for  the  case  of  a  '‘stiff”  tissue. 
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Figure  10:  The  concentric  spheres  problem.  Huge  mesh.  Plot  of  pressure 
along  the  vertical  section  in  dB  for  the  case  of  a  “stiff”  tissue. 
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Figure  11:  The  concentric  spheres  problem.  Big  mesh.  Plot  of  real  part  of 
pressure  on  plane  y  =  0  passing  through  origin. 
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Im(p) 


Figure  12:  The  concentric  spheres  problem.  Big  mesh*  Plot  of  imaginary 
part  of  pressure  on  plane  y  =  0  passing  through  origin. 


16 


Figures  13  and  14  display  the  distribution  of  real  and  imaginary  parts 
of  the  pressure  over  the  =  0  section  obtained  on  the  big  mesh  rescaled  to 
values  from  -.00001  to  .00001  for  the  real  part*  and  from  -.0001  to  .0001  for 
the  imaginary  part. 


Figure  13:  The  concentric  spheres  problem.  Big  mesh.  Plot  of  real  part 
of  pressure  on  plane  y  —  0  passing  through  origin,  in  the  range  -.00001  to 
.00001. 

The  statistics  for  the  largest  mesh  on  which  we  have  succcsfully  run  the 
code  is  as  follows. 


number  of  tets 

=  890144 

number  of  tissue  tets 

=  183872 

number  of  skull  tets 

=  353136 

number  of  cochlea  tets 

=  0 

number  of  air  tets 

=  353136 

number  of  PML  prisms 

=  16816 

total  number  of  d.o.f. 

=  430566 
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Figure  14:  The  concentric  spheres  problem.  Big  mesh.  Plot  of  imaginary 
part  of  pressure  on  plane  y  =  0  passing  through  origin,  in  the  range  -.0001 
to  .0001. 
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3.2.2  Scattering  of  a  plane  wave  on  a  ‘^head^’  with  a  cochlea 

The  analysis  domains  is  divided  into  three  subdomains  as  shown  in  Figure 
15:  cochlea,  skull,  and  air.  The  skull  is  represented  by  a  spherical  shell, 
and  the  cochlea  is  located  inside  the  ear  canal  represented  with  a  cylindrical 
cavity.  Except  cochlea  and  skull,  the  remaining  part  of  region  within  the 
outer  bounding  sphere  is  classified  as  air.  The  cochlea  is  connected  with  the 
skull  by  '"growing”  manually  an  additional  bone  structure  in  between  the 
skull  and  the  cochlea,  i.e.  a  number  of  air  elements  in  the  cavity  is  manually 
reclassified  as  the  skull  elements. 


Figure  15:  The  analysis  domain  of  the  human  hearing  system  with  three 
materials:  cochlea  (yellow),  skull  (green)  and  air  (blue). 

The  head  is  excited  with  the  same  plane  wave  as  in  the  first  example. 
The  purpose  of  this  example  was  to  investigate  the  proportion  of  the  energy 
transferred  to  the  cochlea  directly  through  air  and  indirectly  through  the 
bone.  It  took  several  iterations  to  generate  a  FE  mesh  that  captures  geo¬ 
metrical  details  of  the  coclilea,  provides  a  iiessary  resolution  to  resolve  the 
wavelength  scales  in  the  “brain” ,  and  satifics  mesh  regularity  criteria  for  the 
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FE  computations  (sulEciently  large  jacobians).  The  statistics  of  the  mesh  is 


as  follows* 


total  number  of  elements  =  955234 
number  of  tets  =  878314 
number  of  tissue  tets  =  260800 
number  of  skull  tets  =  304333 
number  of  cochlea  tets  ==  13168 
number  of  air  tets  =  300013 
total  number  of  d,o.f,  =  477961 


(3.3) 


A  cross-section  through  the  mesh  and  details  of  the  meshes  of  the  ^'car 
channel’’  and  the  cochlea  model  are  shown  in  Figures  16,  17  and  18* 


Unfortunately,  we  have  not  managed  to  succesfully  run  the  problem  using 
the  MUMPS  solver  due  to  memory  limitations  {the  interface  with  the  parallel 
solver  has  not  been  updated  yet  to  quadratic  elements). 
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Figure  17:  The  “eochlca  problem".  A  detail  of  the  mesh  in  the  “car  channel” 
around  the  model  of  cochlea. 


Figure  18:  The  “cochlea  problem".  Mesh  for  the  cochlea  model. 
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4  Results  of  Phase  I 


The  following  is  a  short  summary  of  our  accomplishments  in  Phase  I  of  tlic 
project, 

1.  We  have  developed  a  mathematical  formulation  for  the  coupled  prob¬ 
lem  and  the  corresponding  Galerkin  approximation  based  on  hybrid 
tctrahcdral/prismatic  meshes, 

2.  We  have  developed  from  scratch  a  FE  code  to  solve  the  coupled  prob¬ 
lem  using  linear  and  quadratic  tetrahedral  elements,  and  prisms  of 
variable  order  in  the  radial  direction  to  handle  the  PML  truncation 
of  the  computational  domain*  The  code  includes:  mesh  generator, 
element  routines,  an  interface  with  VTK  visualization  package  and  in¬ 
terfaces  with  three  linear  solvers;  an  in-house  frontal  solver,  European 
MUMPS,  and  a  parallel  solver.  The  code  consists  of  over  30,000  lines, 
excluding  the  solvers, 

3.  W''e  have  verified  the  code  by  using  the  method  of  manufactured  solu¬ 
tions  and  comparing  results  obtained  with  different  solvers. 

4*  We  have  succesfully  solved  the  concentric  spheres  problem  for  the  case 
of  a  “stiff  tissue”  and  demonstrated  a  convergent  solution  by  employing 
three  different  meshes  with  an  increasing  number  of  d.O-f. 

5,  We  have  generated  meshes  for  the  cochlea  problem  but  have  not  man¬ 
aged  to  solve  it* 

6.  We  have  generated  meshes  for  the  actual  head  problem  that  have  been 
used  by  the  partners  at  Monopole. 

5  Lessons  learned  from  Phase  I 

Our  initial  plans  were  to  use  our  existing  three-dimensional  /ip-code  based 
on  hexahedral  meshes.  The  code  has  been  in  use  for  many  years,  and  it 
comes  with  the  possibility  of  running  various  adaptive  schemes  to  verify  the 
convergence.  Unfortunately,  generation  of  reasonably  regular  hexahedral 
meshes  for  the  problem  of  interest  with  our  existing  software,  has  turned 
out  to  be  unrealistic.  We  had  to  build  a  completely  new  code  for  tetrahedral 
meshes  from  scratch. 

A  lesson  learned  from  the  concentric  spheres  problem  is  that  we  will  need 
higher  order  elements  to  cope  with  the  large  material  contrast. 
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Experience  with  the  cochlea  example  indicates  the  need  for  a  strict  mon¬ 
itoring  of  mesh  quality.  The  scaling  analysis  indicates  that  the  double  pre¬ 
cision  should  t)e  sufficient  for  solving  the  problem  with  a  direct  multifrontal 
solver.  The  encountered  numerical  singularity  occurred  in  the  region  of 
small  elements  surrounding  the  cochlea  model  and  indicates  a  conditioning 
problem  most  likely  related  to  a  local  loss  of  shape  regularity  of  generated 
elements. 
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Appendices 


A  Formulation  of  the  Coupled  Elasticity /Acoustics 
Problem 

In  this  opening  section,  wc  review  the  derivation  of  the  variational  formu¬ 
lations  for  acoustics,  elasticity  and  then  for  the  ultimate  coupled  elastic¬ 
ity/acoustics  problem. 

A,1  Linear  Acoustics  Equations 

The  classical  linear  acoustics  equations  are  obtained  by  linearizing  the  isen- 
tropic  form  of  the  compressible  Euler  equations  expressed  in  terms  of  den¬ 
sity  p  and  velocity  vector  Vi^  around  the  hydrostatic  equilibrium  position 
p  —  po»  t'i  —  0.  Perturbing  the  solution  around  the  equilibrium  position, 

P  =  Po  +  Vi  =  0  +  5vu 

and  linearizing  the  Euler  equations,  see  c.g.  [8|,  we  obtain  a  system  of  four 
first  order  equations  in  terms  of  unknown  perturbations  of  density  6p  and 
velocity  Svi, 

(  (M,(  +  Po{Svj),j  =  0 

I  +  (M,i  =  0  . 

with  6p  denoting  the  perturbation  in  pressure.  For  the  iscntropic^  flow,  the 
pressure  is  simply  an  algebraic  function  of  density, 

p  ^  pip) 

Linearization  around  the  equilibrium  position  leads  to  the  relation  between 
the  perturbation  in  density  and  the  corresponding  perturbation  in  pressure 

P  =  p(^  +  ^{po)^P 

pa 

^The  entropy  is  assumed  to  be  constant  throughout  the  whole  domain 
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Here  po  is  the  hydrostatic  pressure^  and  the  derivative  ^{pn)  is  interpreted 
a  posteriori  as  the  sound  speed  squared,  and  denoted  by  Consequently, 
the  perturbation  in  pressure  and  density  are  related  by  the  simple  linear 
equation, 

6p  —  (?&p 

It  is  customary  to  express  the  equations  of  linear  acoustics  in  pressure  rather 
than  density.  Dropping  deltas  in  the  notation,  wc  obtain, 

r  =0 

I  PoVi^t  +P,t  =0 

In  this  report,  we  shall  consider  only  time-harmonic  problems.  Assuming 
ansatz, 

p((,  a:)  =  e^^p[x),  Ui((,  x)  =  (x) , 

WC  reduce  the  acoustics  equations  to, 

{c“^iwp  +paVjj  =  0 

poiuivi  +p,i  =  0 


or  in  the  operator  form, 


c  ^iwp  +poV  ■  V  =  0 
Poi^Vi  +  Vp  =  0 


(A.1) 


Eliminating  the  velocity,  wc  obtain  the  Helmholtz  equation  for  the  pressure, 

— Ap  —  k^p  —  0  , 
with  the  wave  number  k  “  uj/c. 

Having  obtained  the  second  order  problem,  wc  can  proceed  now  with 
the  derivation  of  the  weak  formulation,  as  it  is  usually  done  in  most  of  text 
books  on  the  subject*  It  is  a  little  more  iluminating  to  obtain  the  same 
variational  formulation  starting  with  the  first  order  system.  First  of  all,  we 
make  a  clear  choice  in  a  way  we  treat  the  two  equations.  The  equation  of 
continuity  (conservation  of  mass)  is  going  to  be  satisfied  only  in  the  weak 
sense,  i.e.  we  multiply  it  with  a  tet  function  g,  integrate  over  domain  fl 
and  integrate  the  second  term  by  parts  to  obtain, 

~  dx-\-  pQ  j  Vnq  dS  =  0,  (A. 2) 
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Here  Vn  =  VjUj  denotes  the  normal  component  of  the  velocity  on  the  bound¬ 
ary. 

The  second  equation  (conservation  of  momentum)  is  satisfied  in  the 
stivng  sense^  i.e.  pointwise.  Solving  for  the  velocity,  we  get, 


V 


(A.3} 


In  particular,  the  normal  component  of  the  velocity  is  related  to  the  normal 
derivative  of  the  pressure, 


1  dp 
= - 

pottjj  on 

At  this  point  we  introduce  different  boundary  conditions: 


•  a  soft  boundary  F d, 


P^Po 


•  a  hard  bouiidai'y  Tjsf, 

Vn  =  Vo 

•  and  an  impedance  condition  with  a  constant  d  >  0, 

rjn  =  dp  +  Vo 

Multiplying  Equation  A. 2  with  iu;,  substituting  the  boundary  data  into  the 
boundary  term,  and  eliminating  the  velocity  in  the  domain  integral  term, 
using  formula  A, 3,  we  get  the  final  variational  formulation. 

p  =  po  on  Fd 

^  )  pg]  dx-\ribjpod  [  pqdS  =  —  /  voqdS 

Jn\  ^  dr^rurc 

Wg  :  (j  =  0  on  Fd 

(A.4} 

We  have  obtained  the  weak  formulation  without  introducing  the  second  or¬ 
der  problem  at  all!  We  have  a  clear  understanding  which  of  the  starting 
equations  is  understood  in  the  weak,  and  which  in  a  strong  sense.  The 
momentum  equations,  consistently  with  their  pointwise  interpretation,  have 
been  extended  to  the  boundary  to  yield  the  appropriate  boundary  condi¬ 
tions,  We  mention  only  that  all  these  considerations  can  be  made  more 
precise  by  introducing  the  language  of  distributions  and  Sobolev  spaces. 
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A. 2  Linear  Elasticity 

The  time-hai'monic  linear  elasticity  equations  include: 

•  balance  of  momentum, 

—pu)  Ui  —  <JijJ  =  fi 

•  Cauchy  displacement-strain  relation, 

^ij  ^  +  ^j,i) 


•  consititutive  law, 

The  tensor  of  elasticities  satisfies  the  usual  symmetry  assumptions, 

Eijkl  ~  ^ijkl  =  Eijki  =  Ekiij  . 

In  the  case  of  an  isotropic  material, 


Eijki  =  P‘(5ik&ji  +  Sii6jk)  + 

and  the  constitutive  law  reduces  to  the  Hooke’s  law, 

Utilizing  the  Cauchy  geometric  relations,  we  eliminate  the  strain  tensor  and 
represent  the  stresses  directly  in  terms  of  the  displacement  gi-adient, 

=  Eijkiukj  (A.  5) 

or,  for  the  Hooke’s  law, 


(Tij  =  +  Xu^,kSij 


(A,6) 


The  momentum  equations  will  be  satisfied  in  the  weak  sense.  We  multiply 
them  with  a  test  function  Ui,  integrate  over  0  and  integrate  by  parts  to 
obtain, 


/  "  puJ^UiVi)  dx  -  aijUjVi  dS  = 

Jii  Jr 


/  fiVidx, 

Jn 


Vvi 


We  introduce  now  the  boundary  conditions. 


(AJ) 
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•  prescribed  displacements  on  F o 


Ui  -  Ui,D 


•  prescribed  tractions  on 


t{  CijUj  —  Qi 


•  prescribed  impedance  on  Fcy 


ti  +  PijUj  ^  Qi 


We  restrict  ourselves  now  to  Wj  =  0  on  Fij^  substitute  the  boundary  data 
into  the  boundary  term  in  Equation  A. 7,  to  obtain, 


■3?  +  /  9iVi  dS 

Jr ^ ur c 

\fvi  :  Ui  =  0  on  Fd 


The  final  variational  formulation  is  obtained  by  substituting  formula  A. 5  for 
stresses, 


/  9iVi  dS 

J  rpfKjTc 

Vvi  :  t;x  =  0  on  F/j 
(A.8) 


Wc  record  the  final  foraulas  for  the  bilinear  and  linear  forms* 


A,3  Elasticity  Coupled  with  Acoustics 

Let  il  be  a  doniaiu  in  In  the  following  discussion  wc  shall  assume  that 
the  domain  H  is  bounded.  We  assume  that  SI  is  split  into  two  disjoint  parts: 
a  subdomain  occupied  by  a  linear  elastic  medium,  and  a  subdomain  ^Iq 
occupied  by  an  acoustical  fluid*  The  two  subdomains  are  separated  by  an 
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interface  F/.  Neither  the  subdomains  nor  the  interface  need  to  be  connected 
(they  may  consist  of  several  separate  pieces).  The  external  boundary  dQ 
will  be  partitioned  into  Diriciiletj  Neumann  and  Caudiy  parts:  r£),riv,rcr 
respectively.  Each  of  these  boundary  parts  may  consist  of  a  part  belonging 
to  the  boundary  5Sle  of  the  clastic  subdomain,  or  the  boundary  d^la  of  the 
acoustical  subdomain.  Using  a  more  precise  mathematical  language, 
axe  assumed  to  be  opened  and  disjoint  and, 

Similarly,  elastic  and  acoustics  parts  of  the  Dirichlet  boundary:  FDcFDa, 
of  the  Neumann  boundary:  F/zc,  F  A^a,  and  the  Cauchy  boundary:  Tcey  FCa, 
axe  open  submanifolds  of  90  and, 

^  F£^e  U  Foo  U  FiVe  U  Tno  ^  F^e  U  Fca  , 

as  well  as, 

^  Ti  UTp^UT  =  F/  U  Foa  U  F/Va  U  Fca  • 

A  two-dimensional  illustration  of  the  scenario  is  shown  in  Figure  19. 

The  coupled  problem  involves  solving  linear  elasticity  equations  dis¬ 
cussed  in  Section  A. 2  satisfied  in  subdomain  fig  coupled  with  the  equations  of 
linear  acoustics  discussed  in  Section  A.l  and  satisfied  in  subdomain  fla.  The 
unknowns  include  the  components  of  the  displacement  vector  6  lie 

and  the  acoustical  pressure  p{x),x  €  Qa-  The  two  sets  of  equations  are  ac¬ 
companied  by  appropriate  boundary  conditions  and  coupled  by  the  following 
interface  conditions: 

1  dp 

lUJUiTli  =  ViUi  = - TT — 71^,  tj  =  dijUj  =  —pui 

pfiuj  axi 

The  first  equation  above  expresses  the  continuity  of  normal  component  of 
the  velocity;  the  normal  clastic  velocity  has  to  match  the  normal  component 
of  the  acoustical  velocity.  The  second  equation  expresses  the  continuity  of 
stresses:  the  normal  elastic  stress  must  be  equal  to  the  (negative)  pressure, 
whereas  the  tangential  component  of  the  clastic  stress  vector  is  set  to  zero, 
since  the  fluid  docs  not  support  a  shear  stress.  As  usual,  uj  is  the  aiigulai' 
frequency,  i  is  the  imaginary  unit,  pf  stands  for  the  density  of  the  fluid, 
and  Hi  denote  components  of  a  unit  vector  normal  to  interface  F/  which  we 
will  assume  to  be  directed  from  the  elastic  into  the  acoustical  subdomain. 
Multiplying  the  first  interface  condition  by  p/w,  we  get, 

dp 

Pf uj^v^  =  ^  “P”*  (A.  10) 


29 


Figure  19:  Topology  of  a  coupled  problem 


where  Uji  =  uiiij  denotes  the  normal  displacement.  Prom  the  mathematical 
point  of  view,  the  conditions  of  this  type  arc  classified  as  weak  coupling 
conditioTis.  The  word  ^‘weak”  refers  here  to  the  fact  that  the  primary  variable 
for  elasticity  -  the  displacement  vector,  matches  the  secondary  variable  (the 
flux)  for  the  acoustic  problem  -  the  normal  velocity  which  is  related  to 
the  normal  derivative  of  pressure.  Conversely,  the  primary  variable  for  the 
acoustic  problem  -  the  pressure,  defines  the  flux  for  the  elasticity  problem. 
This  “cross-coupling'’  is  very  essential  in  proving  the  wcll-poscdncss  of  the 
problem,  and  stabihty  of  Galerkin  approximatioES. 

On  top  of  the  interface  conditions  we  have  the  usual  boundary  conditions 
for  acoustics, 

•  prescribed  pressure  on 


P  =  PD 
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•  prescribed  normal  velocity  on  F/vaj 


t^n  =  t;o 


•  an  impedance  condition  with  an  impedance  constant  d  >  0  on  Fcn, 

Vn  =  dp  +  VQ 


and  for  the  elasticity, 

•  prescribed  displacements  on  Fdci 

Ui  =  Ui^D 


•  pro45cribcd  tractions  on  Fyvo 


tl  : —  (Jijflj  —  Qi 


•  prescribed  impedance  on  Fce, 

ti  +  iw0ijUj  =  gi 


Wc  proceed  now  with  the  derivation  of  the  variational  formulation*  We 
start  with  the  weak  form  of  the  continuity  equation  for  acoustics, 

X  ~  j  +  Pf  '^9 

and  the  weak  form  of  the  conservation  of  momentum  for  elasticity. 


Ps^  ^t^i)  dx  I 


aijUjVi  dS 


fiVi  dx, 


with  and  fi  denoting  the  density  of  solid  and  body  forces,  respectively* 
Boundary  SSIq  of  the  acoustic  subdomain  is  now  split  into  the  interface  F /v 
and  parts  Foai  F^vai  Fca-  For  the  interface  F/,  wc  use  the  first  interface 
condition  to  replace  the  flux  term  pfVjj  with  iufpfUn,  and  proceed  in  the 
standard  way  with  the  acoustic  boundary  conditions,  to  obtain  the  varia¬ 
tional  statement, 


p  ^  pD  on  T Oa 


f  (^pq  —VpV(i\  dx-\-  [  pfdpqdS+  [  iupf  U^,  q  dS  =  [ 

Jua  \  ^  /  Jfca  Jri  JrNa^rca 


\fq  :  q 


pfvoq  dS, 
0  on  r  Da 
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Similarly,  boundary  dile  of  the  elastic  subdomain  is  split  into  the  interface 
r /v  and  parts  rjves  For  the  interface  P /,  we  use  the  second  interface 
condition  to  replace  the  flux  term  aijUj  with  —pui,  and  use  the  boundary 
conditions  to  obtain  the  variational  statement, 

u  —  Up  on 

<  /  {EijkiUkjVij  -  pa0\vi)  dx  +  w  /  PijUjVi  dS+  pv„  dS=  fiVi  dx  +  /  giVi  dS, 

Jn,  dvct  Jri  Jfic  JvN^'jrcc 

Vu  :  t;  =  0  on  r^e 

Multiplying  the  variational  statement  for  acoustics  by  factor  iu>,  we  get  the 
final  variational  formulation  for  the  coupled  problem  in  the  form, 


'  u  €  ud  +  V,  p  e  pD  +  V, 

<  bee(u,v)  +bae(p,v)  =  le(v),  Vv  6  V  (A.  11) 

^  6ea('^t  “S"  6aa  (p,g)  =  ia{qf),  VgeV 

where: 


•  the  bilinear  and  linear  forms  are  given  by  the  formulas: 

bee(u,v)  =  /  (EijkiUk,iVij  -  dx  +  iui  /  iSijUjVi  dS 

Jne  Jrcc 

pVfi  dS 


baeiP 


',V)=  f 
Jr, 


bea{u,  q)  =  -uJ^Pf  /  Unq  dS 

Jr, 

baaip,  q)  =  /  (VpVq  -  k'^pq)  dx  +  iui  /  pfd  pq  dS 

-/rcn 

ie(v)  =  /  fiVi  dx+  giVi  dS 

!„((/)=  iupf  voqdS, 

'/rjVaUr<7a 


(A.12) 


•  ud  e  J/'{ne)  ■=  is  a  finite  energy  lift  of  displacements  uo 

prescribed  on  Fcei  P~D  €  is  a  finite  cnery  lift  of  pressure  po 

prescribed  on  r£)a, 


•  V  and  V  arc  the  spaces  of  the  test  functions, 

V  =  {v  e  :  u  =  0  on  roj 

V  ={q^  ■  q  =  (ionTDa} 


{A.13) 
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•  fc  =  ujjc  is  the  acoustic  wave  number. 

Coupled  problem  A.  11  is  symmetric  if  and  only  if  diagonal  forms  bee  ajid 
baa  su:c  symmetric  and, 


i>ac(P,«)  =  f>ea(«,P)  ■ 

Thus,  in  order  to  enable  the  symmetry  of  the  formulation^,  we  need  to  rescale 
problem  by,  for  instance,  dividing  the  second  equation  by  factor 

B  Finite  Element  Discretization  and  Implementa¬ 
tion  Details 

Except  for  the  PML  domain,  both  acoustic  and  clastic  domains  arc  dis¬ 
cretized  with  the  simplest  linear  tetrahedra,  i.c.  pressure  p  and  clastic  dis¬ 
placement  components  Ui  are  linear  within  each  element.  This  implies  that 
all  interfaces  including  the  truncating  sphere  arc  approximated  with  plane 
triangular  panels.  The  triangular  mesh  on  the  (approximate)  truncating 
sphere  is  extended  in  the  radial  direction  to  form  two  layers  of  prismatic 
elements.  In  order  to  approximate  well  the  PML  induced  layer,  higher  order 
polynomials  in  the  radial  direction  arc  used,  /?  =  4  in  the  first  layer,  and 
p  —  2  in  the  second  layer.  This  is  in  accordance  with  our  experience  of 
resolving  PML  induced  boundary  layers  with  /ip-adaptive  elements,  see  [9] 
for  examples, 

B,1  Generation  of  tetrahedral  meshes 

Wc  choose  an  octree-based  isocontouring  method  [16]  to  extract  interior  and 
exterior  tetrahedral  meshes  for  the  acoustic  and  elastic  domains.  First  wc 
define  the  analysis  domains  and  construct  a  signed  distance  map.  Then 
a  top-down  octree  subdivision  coupled  with  the  dual  contouring  method 
is  used  to  rapidly  extract  adaptive  3D  finite  element  meshes  with  correct 
topology  from  the  signed  distance  map.  Finally,  the  edge  contraction  and 
smoothing  methods  arc  used  to  improve  the  mesh  quality.  This  octree- 
based  technique  extends  the  dual  contouring  method  to  crack- free  interval 
volume  3D  meshing  with  feature  sensitive  adaptation.  Compared  to  other 
tetrahedral  extraction  methods  from  imaging  data,  this  method  generates 
adaptive  and  quality  3D  meshes  without  introducing  any  hanging  nodes. 

At  the  end,  the  following  files  are  created: 

^This  is  essential,  among  other  reasons,  from  the  point  of  view  of  using  a  direct  solver. 
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•  sphere-file  -  contains  a  list  of  vertices  and  triangles  on  the  truncating 
sphere, 

•  tissue^skulLfile  -  contains  a  list  of  vertices  and  triangles  on  the  tis- 
sne/skull  interface, 

•  air^skulLfile  -  contains  a  list  of  vertices  and  triangles  on  the  air /skull 
interface, 

•  C£>cfe/ea-,azr^^i€  -  contains  a  list  of  vertices  and  triangles  on  the  cochlea/air 
interface, 

•  cochlea-file  -  contains  a  list  of  vertices  and  tetrahedra  within  the  cochlea, 

•  air-file  -  contains  a  list  of  vertices  and  tetrahedra  within  air, 

•  skull-file  -  contains  a  list  of  vertices  and  tetrahedra  within  skull, 

•  tissue-file  -  contains  a  list  of  vertices  and  tetrahedra  within  tissue 

The  meshes  arc  fully  compatible,  i*e*  for  instance  all  vertices  for  the  skull 
tetrahedra,  that  are  located  on  the  skull/air  interface,  coicide  with  vertices 
listed  in  aiv-skulLfile. 

B.2  Data  structure  and  element  computations 

An  existing  data  structure  for  higher  order  hcxahcdral  elements,  see  [3,  4], 
has  been  extended  to  the  case  of  tetrahedral  and  prismatic  elements.  The 
data  structure  aiTays  are  initiated  with  a  relevant  information  on  nodal 
connectivities,  and  element  neighbors  necessary  for  element  computations. 

Element  matrices  corresponding  to  bilinear  forms  are  integrated  using 
standard  Gaussian  quadrature  for  tetrahedra  (volume  integrals)  and  trian¬ 
gles  {interface  terms). 

B,3  Solvers 

Three  linear  solvers  are  used  in  this  project.  The  first  one  is  a  serial  frontal 
solver  developed  at  ICES^,  the  second  one  is  the  European  MUMPS,  sec 
[10],  and  the  third  one  is  a  new,  parallel  solver  enabling  solution  of  large 
systems  of  equations  with  several  milions  of  unknowns,  developed  especially 
for  this  project.  Preliminary  details  on  the  solver  are  given  in  Section  C,  The 

^The  solver  was  developed  by  Dr.  Eric  Beeker,  a  professor  in  the  ASE/EM  Dept,  and 
a  long  time  member  of  TICOM,  next  TICAM  and  now  ICES 
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interfaces  witli  the  first  two  solvers  have  been  implemented  for  the  purpose 
of  verification, 

B.4  Graphics 

To  visualize  our  models,  meshes  and  solutions,  we  have  implemented  a  simple 
interface  to  the  Visualization  Toolkit  {VTK)  [11],  a  collection  of  C++  classes 
that  implement  a  wide  range  of  visualization  algorithms.  The  central  data 
structure  for  this  project  is  the  vtkUnstructuredGrid,  which  represents 
volumetric  data  as  a  collection  of  points  with  corresponding  scalar  values 
(in  this  case,  the  real  and  imaginary  part  of  the  pressure),  connected  by 
cells  of  arbitrary  type  and  dimension  (i.e,  lines,  triangles,  quads,  tetrahedra, 
prisms,  etc,).  This  data  is  then  plugged  into  a  variety  of  filters  that  allow 
us  to,  for  example,  “slice’"  through  the  dataset  with  a  plane  to  sec  the 
pressure  in  the  interior,  extract  colored  isocontours  or  isosurfaces  of  the 
pressure,  or  generate  an  animation  of  the  time-dependent  pressure  P(x,  i)  — 
Rc 

C  Parallel  Linear  Solver 

The  nested-dissections  parallel  multi-frontal  solver  is  utilized  to  solve  the 
problem  over  large  meshes.  The  frontal  solver  is  an  extension  of  the  Gaus¬ 
sian  elimination,  where  assembly  and  elimination  arc  performed  together  on 
the  so-called  frontal  sub-matrix  of  the  global  matrix  [7],  The  multi- frontal 
solver  utilizes  domain  decomposition  pattern  to  work  with  multiple  frontal 
matrices. 

An  example  of  the  computational  domain  related  to  the  head  problem, 
partitioned  into  3  sub-domains  is  presented  in  Fig.  20.  Local  orderings 
of  degrees  of  freedom  (d.o.L)  are  computed  on  each  sub-domain,  where 
internal  d.o.f.  are  numbered  first,  and  interface  d.o.f.  arc  put  at  the  end. 
Local  matrices  contain  internal  d*o.L  interactions  part  A;,  interface  d.o.f. 
interaction  part  and  the  interface- internal  d.o.f.  interaction  parts  Ri, 


'  .4; 

Bi  ' 

Xi 

■  bi  ' 

The  global  matrix  is  a  sum  of  local  matrices 

d 

A  =  (C-2) 

i=l 
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Figure  20:  Computational  mesh  partitioned  into  3  sub-domains. 

The  global  numbering  of  d*o.L  is  created  in  the  following  way: 

•  All  internal  d.o.f.  from  1-st  sub-domain, 

•  All  internal  d.od.  from  2-nd  sub-domain, 

•  .*♦ 

♦  All  internal  d.o.fr  from  d-th  sub-domain, 

•  All  d.o.fr  from  the  global  interface. 

The  maps  Pi  arc  constructed  to  transfer  ordering  from  local  to  global  Going 
into  details,  let  n  =  denotrs  total  number  of  d.oi., 

^[nternai  number  of  interior  d.o.f.  of  sub-domain  L 

^interface  ”  number  of  interface  d.o.f.  of  sub-domain  i. 

Hinterface  =  global  number  of  interface  d.o.l. 

We  define  the  P,  maps  in  the  following  way: 

^  ("internal  +  "ilrface  "inlernnl  +  "ilrfnce)  ^  M  [n  ^  n)  (C,3) 

^  ("flrface  ><  "Lrnnl)  ^  ("llrface  ^  "i?terfnce)  (^-4) 
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Figure  21:  Forward  elimination  is  stopped  before  processing  fake  elements* 
Sub-domain  internal  nodes  are  eliminated* 


Pbb  -  ^  (^Lterface  ^  ^LLrfoce)  ^  (’T-interface  ^  Tljnterface) 


Since 


PiAiPl  = 


0 


Ai  PibBiP^ 

PlnCiPZ  Pb,A^,Pl 

the  structure  of  the  global  stiffness  matrix  is 


‘  Xi  “ 

A2 


Aa 


PibBiPfb 

PibB2P^ 


(C.5) 

(C.6) 


{C,7) 


"  XI  ^ 

X4 

.  ^3  _ 

nbf>i  Pfb 

P2t>hP£ 

P.iijl>dP]l 

EL  Pitb.Pl 

(C.8) 


In  the  3D  code  we  utilize  a  serial  version  of  the  MUMPS  solver  [10] 
over  each  sub-domain.  Single  processor  versions  of  MUMPS  arc  executed 
on  each  processor  for  each  sub-domain*  The  MUMPS  is  executed  to  provide 
the  Schur  complement  of  the  local  sub-domain  matrix.  The  internal  sub- 
domain  d.o.f*  are  eliminated  with  respect  to  the  interface  sub-domain  d*o.f. 
The  Schur  complement  of  the  local  system  (C.l)  for  sub-domain  i  is 

=  bi  (C*9) 
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Figure  22:  Global  interface  problem  can  be  aggregated  by  summing  up  local 
Schur  complements* 


Ai  =  Ai~CiA-^Bf  (C.iO) 

=  (C-U) 

A  practical  way  of  obtaining  the  Schur  complement  is  to  order  d*oi.  in 
such  a  way  that  interface  d*oi,  are  numbered  at  the  end,  and  than  to  execute 
forward  elimination  for  the  internal  d.o.f.  [13]*  The  forward  elimination  is 
stopped  before  processing  interface  d*o*f.  After  eliminating  local  internal 
nodes,  the  local  matrices  look  like  it  is  presented  in  Fig.  21  and  in  (C.12). 


■  Ui  b;  ■ 

'  Xi  ' 

1 

.  0  . 

Once  local  Schur  complements  are  computed,  we  can  formulate  and  solve 
the  global  interface  problem.  The  Schur  complement  matrices  can  be  sent 
to  separate  process,  the  global  interface  problem  matrix  can  be  obtained  by 
summing  up  the  local  interface  matrices,  sec  Fig.  22. 


II 

(C.13) 

II 

.Sj 

(C.14) 

t=i 

(C.15) 

i-l 


The  global  interface  problem  can  be  solved,  by  utilizing  sequential  version 
of  MUMPS,  see  Fig.  23 

Ax  =  b^Ux  =  L-^b  (C.16) 
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Figure  23:  Global  interface  problem  is  solved* 


Figure  24:  Solution  of  the  global  interface  problem  can  be  broadcasted  into 
sub-domains. 


and  the  solution  can  be  broadcast  into  all  sub^domainsT  see  Fig.  24.  The 
part  of  the  solution  corresponding  to  local  interface  node  is  substituted  to 
the  right-hand-side,  the  identity  matrix  is  put  into  the  right-bottom  part  of 
the  local  matrix 


'  Ui  B*  - 

r  1 

Xi 

1 - 

I 

i _ 
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Xi 

.  . 
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bi 


(C.17) 


and  the  backward  substitution  is  executed  on  sub-domains* 

Summing  up  the  presented  scheme,  to  solve  the  problem  distributed  into 
sub-domains,  wc  need  to  formulate  and  solve  the  interface  problem  first. 
For  large  3D  problems,  or  for  multiple  sub-domains,  the  interface  prob¬ 
lem  is  large,  its  matrix  is  dense.  This  motivated  us  to  utilize  the  idea 
of  nested  dissections  to  reduce  computational  cost  of  the  interface  problem 
solution.  The  idea  of  the  nested  dissection  solver  is  to  utilize  the  Sdmr  com¬ 
plement  pattern  recursively.  This  can  be  described  in  the  following  steps, 
compare  Figure  25. 


•  In  the  first  step,  local  internal  sub-domain  d,o,f.  arc  eliminated  with 
respect  to  local  interface  sub-domain  d,o.b,  in  the  same  way  as  it  is 
described  above* 


•  Processors  arc  joined  into  pairs,  partial  interface  matrices  arc  aggre- 
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Figure  25:  Paatial  Schur  complements  (1-2-3)  followed  by  backward  substi¬ 
tutions  (4-5-6)  executed  over  the  nested  dissections  scheme. 


gated  for  each  pair  of  processors.  Common  interface  d.o.f.  are  elim¬ 
inated  with  respect  to  external  interface  d.o.b  within  each  pair  of 
processors.  Notice  that  this  procedure  requires  to  build  new  inter¬ 
face  node  numbering  within  each  pair  of  processors.  This  is  the  main 
technical  difficulty  in  implementing  the  nested  dissection  scheme. 

•  The  procedure  from  the  previous  step  is  repeated  recursively  as  long 
as  there  are  common  pieces  of  interface  to  be  eliminated. 

•  In  the  last  step  of  the  elimination,  there  is  only  one  piece  of  aggregated 
interface  d.o.f.  matrix,  the  common  interface  part.  This  interface 
problem  matrix  is  much  smaller,  since  it  corresponds  to  common  part 
of  the  interface  between  two  groups  of  processors.  In  general,  when 
the  domain  is  not  cylindrical,  this  interface  part  is  associated  with 
the  cross-scction  of  the  domain.  This  interface  problem  can  be  solved 
now,  since  it  is  much  cheaper  then  the  global  interface  problem  in  the 
previous  scheme. 

•  The  solution  of  the  common  interface  part  is  broadcast  back  into  two 
group  of  processors,  local  solution  extracted  by  using  built  maps  arc 
substituted  into  Schur  complement  matrices  from  current  nested  dis¬ 
sections  step,  and  the  backward  substitution  is  executed  to  get  next 
contribution  to  the  global  interface  problem  solution, 

•  The  procedure  is  repeated  until  we  end  up  with  the  complete  solution 
of  the  global  interface  problem. 
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Figure  26:  Real  part  of  the  solution  on  the  simplified  model,  obtained  by 
parallel  execution  on  5  processors. 


Figure  27;  Imaginary  part  of  the  solution  on  the  simplified  model,  obtained 
by  parallel  execution  on  5  processors. 

•  In  the  last  step,  the  backward  substitution  is  executed  over  sub-domains 
and  the  problem  is  finally  solved. 

The  presented  strategy  can  be  easily  generalized  to  the  number  of  processor 
not  equal  to  the  power  of  2.  The  exemplary  solution  on  the  simplified  model 
of  human  head  build  by  concentric  spheres  is  presented  in  Figures  26  and 
27.  Because  of  material  data  assumed,  the  maximum  solution  values  arc 
obtained  on  the  skull,  compare  28. 

The  current  version  of  the  nested  dissection  parallel  multi- frontal  solver 
attains  relative  efficiency  up  to  60  percent. 
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Figure  28:  Solution  on  the  central  sub-domain,  top  view. 
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